Clustering and fluidization in a one-dimensional granular system: molecular dynamics 

and direct-simulation Monte Carlo method 
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We study a 1-D granular gas of point-like particles not 
subject to gravity between two walls at temperatures Ti e f t 
and Tright- The system exhibits two distinct regimes, depend- 
ing on the normalized temperature difference A = (T r i g ht — 
Tieft) /(Tright + Tieft): one completely fluidized and one in 
which a cluster coexists with the fluidized gas. When A is 
above a certain threshold, cluster formation is fully inhib- 
ited, obtaining a completely fluidized state. The mechanism 
that produces these two phases is explained. In the fluidized 
state the velocity distribution function exhibits peculiar non- 
Gaussian features. For this state, comparison between inte- 
gration of the Boltzmann equation using the direct-simulation 
Monte Carlo method and results stemming from microscopic 
Newtonian molecular dynamics gives good coincidence, estab- 
lishing that the non-Gaussian features observed do not arise 
from the onset of correlations. 



I. INTRODUCTION 

Granular systems have been extensively studied due 
both to the theoretical challenges they present (for a re- 
cent review see [JTJj ) and to the applications of industrial 
importance that spring from the rich phenomena they 
exhibit (see (^|| and references therein). These systems 
are characterized by an energy loss in collisions. This 
loss is at the base of many interesting phenomena, such 
as inelastic collapse where the particles collide in- 

finitely often in finite time, and clustering (for a sample 
of theoretical, simulational, and experimental approaches 
see ^-^] ) . Different methods for keeping the system from 
collapsing have been devised, such as subjecting the par- 
ticles to Brownian forces |lO|Jll| ] , and forcing through the 
boundaries by putting the system in a box with one or 
more thermal- like walls (see for example ]T^dT^|). This 
work focuses on the latter method. 

Being one of the simplest types of forcing, several au- 
thors [^2 13 15|- p7| , ^9[ have studied a one-dimensional 
system in a box with one or two heated (stochastic) 
walls. Of these, |l2|Jl3|Jl^Jl^| study cluster formation, 
although (ljjl^l are not strictly one-dimensional. 

This article studies a quasielastic one-dimensional sys- 
tem not subject to gravity between two thermalizing 



* Electr onic address: jpasinik@cec.uchile.c l 
^URL: http:/ /www.cec.uchile.cl/cinetica/ 



walls. We focus on two control parameters: the total 
inelasticity parameter qN = N(l — r)/2, where N is 
the number of particles and r is the restitution coeffi- 
cient, and the externally imposed temperature gradient. 
The parameter qN has been shown to be relevant for 
the quasielastic system , , 19 - 21 1. By varying these 
parameters we determine the region in parameter space 
where clustering is fully inhibited, obtaining a fluidized 
state. We present a singular feature of the distribution 
function for the clustering regime, and then study how 
this feature is modified for the fluidized state. 

In fl^| the authors study a one-dimensional system 
of point-like particles between an elastic and a heated 
wall. They emphasize that a cluster inevitably forms 
away from the heated wall, regardless of how elastic the 
system is (as long as it is not perfectly elastic) . They also 
study the same system, but with both walls expelling the 
particles with a fixed velocity. In this case they find that 
the cluster forms away from the walls and roams slowly 
about the system, with two groups of fast particles con- 
necting the cluster with the "heated" walls. 

In |13| the same system is studied for different types 
of boundary conditions at the heated wall. The stochas- 
tic boundary condition studied has the form of a power 
of the velocity times the "thermal" condition (the one 
that produces a Maxwell-Boltzmann distribution in the 
elastic case). The authors show that when the power 
that multiplies the thermal condition is positive the test- 
particle equation (derived from the Boltzmann equation) 
has a steady-state solution. Thus the thermal case does 
not have a steady state and develops a cluster away from 
the heated wall. The mechanism for the growth of the 
cluster is explained and verified numerically. 

In a similar system is studied: a long thin pipe 
of inelastic hard disks with heated walls (at the same 
temperature) at the ends of the pipe and periodic side 
walls. The pipe is thin enough for the particle order to 
be preserved. The probability distribution for the dis- 
tance between the central particles is studied. This dis- 
tribution gives a markedly denser system near the center 
than in the elastic case, although the limit to the elastic 
case is smooth, unlike the strictly one-dimensional case 
of |^2] , ^3| . In jll| the same author studies the velocity 
correlations that this system develops as inelasticity is in- 
creased, showing that a consistent description must take 
these correlations into account. 

In this paper we revisit the one-dimensional system 
of N point-like particles interacting via collisions that 
conserve momentum but dissipate kinetic energy. To fix 
notation, the particle velocities after a collision are given 



1 



by 



"i = QCi + (1 - q)c 2 , 4 = (1 - q)c x + qc 2 , 



(1) 



where Ci is the velocity of particle i before a collision, and 
q = (1 — r)/2, being r the restitution coefficient. For the 
elastic case (r = 1) the particles simply exchange veloci- 
ties. Since the particles are point-like, the system is then 
indistinguishable from a system in which the particles do 
not interact. 

This one-dimensional system is interesting because dis- 
sipation is the first order correction to a free gas. Besides, 
results for the one-dimensional system have been found to 
have unexpected relevance for higher-dimensional prob- 
lems. For example, in two-dimensions the particles in- 
volved in inelastic collapse lie roughly on a line ||. Also, 
the dissipation-induced temperature gradients calculated 
m @ for the one-dimensional case inspired the authors 
to look for dissipation-induced Rayleigh-Benard-like con- 
vection for a two-dimensional system without an exter- 
nally imposed temperature gradient [T^ ] . 

For a system with one thermal wall and open on the 
other side, under the influence of gravity, the quasielastic 
system may be kept fluidized [17 19 1: any cluster that 
starts to form is forced against the thermal walk where it 
evaporates. In jl7| the test-particle equation |]l2| , |l3| , pot — 
which is the 1-D Boltzmann equation where the limit 
N — > oo is taken, but keeping qN fixed — is successfully 
applied, with close matching of theory and simulations 
even at the level of the distribution function. 

The one-dimensional system under study is left to 
evolve between two thermal walls at temperatures Ti e f t 
and T r i g ht, with Xi e f t < T r ight- We define the parameter 



T, 



right 



;ht 



Tlcft 



(2) 



to quantify how far from the symmetrical case is the 
system. Under the same conditions, an elastic sys- 
tem has a perfectly bimodal velocity distribution with 
global homogeneous temperature equal to ^/Ti c f t T r ight- 
For the sake of comparison, we simulate systems with 
^/T^cft Tright = 1- Here A = represents a symmetri- 
cal setting, while A = 1 represents an infinitely strong 
temperature gradient. 

As in jlj] , for A = a cluster unavoidably forms away 
from the thermal walls. After forming, the cluster per- 
forms an apparently random walk about the center, grow- 
ing in size (and therefore mass) while the rest of the sys- 
tem grows more rarefied. With the decrease in density 
of the surrounding gas and the increased inertia of the 
cluster, an eventual collision with one of the walls is to be 
expected. When this happens part of the cluster evap- 
orates, and what is left of it is expelled from the wall 
(see Fig. [l]), thus restarting the growth process. Thus 
not only is the system highly clumped, but also in a non- 
steady state. Nevertheless, the gas that is far from the 
random-walk zone has a well-defined time average for the 
distribution function, as is seen in Fig. |[ A noteworthy 




8x10 



FIG. 1. Subfigure (a) shows the time evolution of the num- 
ber of particles inside the cluster for a system of N = 1000 
particles between two walls at the same temperature for 
qN — 0.01. Note the plateau after each fall. Subfigure (b) 
shows the trajectory of the cluster for the same time inter- 
val. The walls are placed at x — and 1. The temperature 
at the walls is unity, and thus a unit of time measures how 
long it takes for a particle with the thermal speed to cross the 
system. 



feature of this distribution is that it exhibits apparently 
singular behavior for slow velocities. 

Setting A the symmetry of the system is bro- 
ken. For A <C 1 the cluster performs a slightly asym- 
metric random walk, spending more time near the colder 
wall, and therefore colliding more often with it. Thus 
the cluster cannot grow as much as it did in the sym- 
metric case before colliding with a wall. By increasing 
A the cluster- wall collision frequency grows, even ob- 
taining short "windows" in which the cluster completely 
evaporates. By further increasing A these windows grow 
larger until a point is reached where no cluster forms. In 
this fashion a totally fluidized state is achieved, which 
may be tractable with the dissipative Boltzmann equa- 
tion. However, the distribution function obtained from 
the molecular dynamics (MD) simulations exhibits a pe- 
culiar non-Gaussian feature for slow velocities. This fea- 
ture is a smoothed version of the apparent singularity 
of the symmetric case. To discern whether this feature 
is due to correlations or is present before they settle in, 
we compared Newtonian molecular dynamics results with 
those obtained through direct-simulation Monte Carlo 
(DSMC), which neglects correlations. The results agree 
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FIG. 2. Subfigure (a) shows the velocity distribution func- 
tion at the left wall for the system referred to in Fig. jjj. The 
distribution for particles reaching the wall (c < 0) exhibits a 
sharp peak. Subfigure (b) shows the same distribution multi- 
plied by yjcf to show that the peak behaves like | c | 1 / 2 _ The 
microscopic velocity c is measured in units of the thermal 
speed. 



very well, except when the system approaches the clus- 
tering regime. 




qN 

FIG. 3. Limiting values of A = (r ri ght - T le{t ) / (T right + T lcit ) 
as a function of qN for N = 1000. y/Ti eftTright = 1 through- 
out. The solid curve shows the lowest possible value A can 
take without detecting clusters. The dashed curve shows the 
largest possible value A can take with easy cluster detection. 
The small region between the curves represents a zone where 
clusters appear erratically. 



are uncommonly near by three orders of magnitude. We 
discard chains of length three or less, since they may be 
random encounters. Measuring the total length of the 
cluster we have found that on average it is of the order 
of 10~ 5 ; thus the choice of 10~~ 5 as link-link distance is 
much larger than the true distance between them. 

As already stated, the boundary conditions are such 
that the (homogeneous) temperature of the correspond- 
ing elastic system (equal to *\/^icft bright) is one. 



III. RESULTS 



II. SIMULATION METHOD 

We simulate the system through event-driven molec- 
ular dynamics 22 and through direct-simulation Monte 
Carlo [p3| . The direct-simulation Monte Carlo proce- 
dures use the null-collision technique p3] where, over- 
estimating the collision frequency (using the maximum 
relative velocity within a cell), the number of collisions 
to be attempted is calculated through a Poisson process. 
In the next step the collisions are attempted, choosing at 
random two particles within the cell, and making them 
collide with a probability proportional to their relative 
velocity. Most of the molecular dynamics and DSMC 
simulations were done with N = 1000 particles. 

In the MD simulations we detect clusters using a geo- 
metric criterion: we consider chains of particles that are 
nearer than a critical distance (in our case 10 -6 and 10~ 5 , 
to be certain that the conclusions are independent of the 
choice). The system length is one, and with a thousand 
particles the mean distance between neighbors for a ho- 
mogeneous system is 10~ 3 . Thus we detect particles that 



A. Clustering regime 

Figure [j] shows the non-steady state of a granular 
system between two walls at the same temperature for 
qN = 0.01. A cluster forms away from the walls, per- 
forming a random walk of varying amplitude. When 
the cluster reaches a wall, part of it evaporates, and the 
growth process begins anew. 

As is usual for the quasielastic case, we relabel the par- 
ticles when they collide. This enables us to visualize this 
system as a group of barely interacting particles passing 
through each other. 

The picture for cluster evolution, as explained before, 
is the following: the cluster grows because the slowest 
particles, due to the asymmetry of the distribution func- 
tion, drift towards the cluster jl|] . As it grows, the den- 
sity of the gas surrounding it decreases, with the conse- 
quent saturation in growth. Thus we have a "Brownian 
particle" of increasing mass moving in an increasingly 
rarefied medium. This "particle" will be increasingly 
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FIG. 4. Number density, temperature, and heat flux pro- 
files for N = 1000, qN = 0.1, and Ti cft = 0.5 (A = 0.6). 
The solid line represents results from a molecular dynam- 
ics simulation, while the dots represent results from a Monte 
Carlo simulation. The density and temperature are related 
by p = nT and, since momentum is conserved and the system 
is stationary, the pressure is constant throughout the system. 
The density and temperature profiles are almost symmetrical 
because the normalized density is very close to unity. 



less affected by the surrounding medium, until it can no 
longer be kept away from the walls. 

The cluster moves several orders of magnitude slower 
than the thermal speed (three orders of magnitude in 
Fig. [I]) . Upon reaching a wall, the front liners strike the 
wall and are expelled by it much faster than the other 
cluster members. These particles pass through the clus- 
ter, transferring momentum to it, as described in [fL2[ . 
Thus these fast particles push the cluster away from the 
wall, where it can absorb particles again. The fast parti- 
cles, however, no longer belong to the cluster. 

Since the slowest particles in the gas are the ones that 
will be absorbed by the cluster it is the number of 
slow particles in the gas that will determine the clus- 
ter's growth rate. After the cluster strikes a wall, the 
expelled particles will be fast particles, and they will not 
contribute to the growth of the cluster during the time it 
takes for the gas to cool down again: the only particles 
available for absorption are the ones that were available 
before the cluster-wall collision. This explains why the 
cluster keeps growing at approximately the same rate it 
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FIG. 5. Velocity distribution function at the cold boundary 
(left wall) for the system referred to in Fig. ^. The solid 
line represents results from a molecular dynamics simulation, 
while the dashed line represents results from a Monte Carlo 
simulation. 



did before the collision. After the gas has cooled down, 
the growth rate returns to its normal value. This is the 
end of the plateau seen in Fig. |l] after each cluster-wall 
collision. 

To quantify the evaporation process we proceed as 
in as soon as the first particle belonging to the 

cluster reaches a wall, it is expelled with a speed much 
higher than the cluster velocity; thus we may consider 
the ideal situation of a cluster of M particles at rest be- 
ing stricken by a fast particle with velocity v (in this case 
v « 1). After colliding with the first particle in the clus- 
ter, the new fast particle's speed will be (1 — q). Thus, 
after traversing the cluster, the fast particle's velocity 
will be (1 — q) N . Since momentum is conserved in colli- 
sions, the center of mass of the cluster will have acquired 
a speed of 
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(3) 



Considering the case N ^> 1 with fixed qN, as in fllTf! , we 
may simplify this expression to 



-qN 



VCM 



N 



(4) 



By further considering the case qN <C 1 we get vqm ~ 9- 
In this limit, if the cluster reaches the wall with veloc- 
ity vq and is expelled from it with velocity v\ , the number 
m of particles evaporated will satisfy \v\ — Uq| « mq. For 
the situation shown in Fig. [l], Vq and V\ are typically of 
the order of 10~ 3 and q = 10~ 5 , hence the number of 
particles evaporated will be of the order of 10 2 . 

Even if the system is in a non-steady state, the gas 
at the walls (far from the random-walk zone) has a well- 
defined time average for the distribution function. The 
distribution function at the left wall is shown in Fig. |^. 
There is an apparent singularity for slow velocities. The 
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distribution is asymmetric as it should, since the particles 
leaving the wall (c > 0) follow a Gaussian distribution. 
Figure || also shows the distribution multiplied by y/\c~\- 
Since the limit of \/[c[/ for c — > CP is finite and nonzero, 
we conclude that the distribution function exhibits a sin- 
gularity that behaves like | c | 1 / 2 for slow velocities. 

As shown in p3| , when A = (the symmetric case) 
the distribution shown in Fig. || is not a solution of the 
steady-state test-particle equation: 



c d x f(x, c, t) = qN d c [f(x, c, t)M(x, c, t)] 



(5) 



where 

M(x,c, t) 



f(x,c\t)(c — c)\c — c'\dc. (6) 



To establish this, let us study the behavior for small c 
of a solution of this equation. Assume that, for small c, 
f(x,c « 0) Rs fo(x)c~ a , with a < 1 in order to have a 
finite density in the vicinity of x. Furthermore, assume 
that M(x,c «0)« Mo(x)c l3 . Inserting this behavior in 
Eq. (|J) we obtain 



d x f « qNM f (P 



„/3-2 



(7) 



Thus, in order to keep /o (the amplitude of the singular- 
ity) finite, we must have either (3 = a or j3 > 2. Inte- 
grating the distribution of Fig. ^ we obtain (3—1. Since 
(3 < 2, we must have a = 1. But this corresponds to a 
nonintegrable distribution, and therefore the distribution 
cannot be steady. 



B. Inhibition of cluster formation 

Figure || shows the regions in (A, qN)-spa.ce where 
clustering is inhibited for N = 1000. As is to be ex- 
pected, as the inelasticity increases, a stronger tempera- 
ture gradient is necessary to inhibit cluster formation. 

To discern whether the non-Gaussian features of the 
velocity distribution function are derived from correla- 
tions in the system we compared results from MD simula- 
tions (full Newtonian dynamics) with results from DSMC 
simulations (no velocity correlations assumed). Figures |] 
and |^ show this comparison for a case far from the cluster- 
ing threshold (A = 0.6 and qN = 0.1). The temperature 
of the left and right walls are chosen so that the global 
temperature for the elastic case (equal to y / TTJftTright) is 
one. The curves match almost exactly. 

At the level of the distribution function, the results 
also match closely. The peculiar non-Gaussian feature of 
the distribution function is clearly seen in Fig. |. There 
is some slight mismatch near the peak. 

For the fluidized case, since momentum is conserved 
and the system is stationary, the pressure is constant 
throughout the system. The number density and the 
granular temperature calculated here are related by p — 
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FIG. 6. Number density, temperature, and heat flux pro- 
files for N = 1000, qN = 0.1, and Ti cft = 0.84 (A = 0.173). 
The solid line represents results from a molecular dynam- 
ics simulation, while the dots represent results from a Monte 
Carlo simulation. 



nT (T in energy units). Thus when the normalized den- 
sity n/N varies little throughout the system (n/N en 
1 + e(x)), the normalized temperature is 



1 



T p/n tiq/N 
To p/n Q n/N 1 + e(x) 



l-e(x), (8) 



thus obtaining the nearly symmetric profiles seen in 
Figs. | and E 

Figures |[ 0, and || compare the MD and DSMC results 
for cases near cluster formation. The non-Gaussian fea- 
ture of the distribution function shows a systematic devi- 
ation for DSMC simulations: there is overpopulation for 
slow velocities. This is explained by considering that the 
DSMC method, like the Boltzmann equation, neglects 
correlations. When the system approaches the clustering 
regime, increased dissipation induces correlations which 
tend to make the particles collide less 16 1. In DSMC 



these correlations are neglected, with the corresponding 
systematic overestimation in the collision frequency. This 
overestimation results in a lower temperature of the sys- 
tem about the density peak. 
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FIG. 7. Number density, temperature, and heat flux pro- 
files for N = 1000, qN = 0.5, and Ti cft = 0.01 (A = 0.9998). 
The solid line represents results from a molecular dynam- 
ics simulation, while the dots represent results from a Monte 
Carlo simulation. 



IV. CONCLUSIONS 

We have shown that a system not subject to gravity 
between thermal walls unavoidably reaches a non-steady 
state when the walls are at the same temperature. A clus- 
ter forms in the bulk, slowly roaming about the system 
while absorbing particles. As it grows, the amplitude of 
the random walk increases, until at last the surrounding 
gas cannot keep the cluster away from the walls. When 
the cluster reaches a wall, a part of it is ejected by the 
wall through the cluster (relabeling the particles on colli- 
sions), effectively pushing the cluster away from the wall, 
and leaving it to grow again. 

Most of the time the cluster is far from the walls. Thus 
measuring the distribution function at a wall is measur- 
ing the distribution function of the gas that surrounds 
the cluster. This distribution function has a well-defined 
time average, and exhibits apparently singular behavior 
for slow particles, diverging like cl" 1 / 2 . 

Imposing an external temperature gradient forces the 
cluster against the colder wall, inhibiting its growth. In- 
creasing the temperature difference leads to a system in 
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FIG. 8. Velocity distribution functions at the cold bound- 
ary (left wall). Subfigure (a) corresponds to the system re- 
ferred to in Fig. [| while subfigure (b) corresponds to the 
system of Fig. m The solid line represents results from a 
molecular dynamics (MD) simulation, while the dashed line 
represents results from a Monte Carlo (DSMC) simulation. 
The MD results have been rescaled so that the area under 
the MD and DSMC curves is the same. There is a systematic 
overpopulation of slow particles in the distributions obtained 
from Monte Carlo simulations. 



which the cluster never forms: the system is completely 
fluidized. The distribution function of the gas exhibits 
peculiar non-Gaussian features: a smooth version of the 
aforementioned singularity. Therefore, any attempt at 
solving the Boltzmann equation through moment meth- 
ods must consider this feature in the initial ansatz, as is 
done for the problem of an infinitely strong shock wave 
in p5| and pq| . In fact, a solution for this problem was 
attempted using the four moment method of |^7| . As 
mentioned in 28 1, the fourth balance equation could not 
be freely chosen when the boundary conditions were sym- 
metric: some choices gave undefined results. As is nat- 
ural, by not including the non-Gaussian feature in the 
ansatz for this calculation we obtained absurd results, 
such as higher temperature in the middle of the system 
than near the walls. 

We compared molecular dynamics with direct- 
simulation Monte Carlo. Agreement between these two 
methods shows that the non-Gaussian feature of the dis- 
tribution function may be predicted by the dissipative 
Boltzmann equation. As the system approaches clus- 



G 



ter formation, correlations settle in. These correlations 
reduce the collision frequency among particles. DSMC 
neglects these correlations, and thus overestimates the 
number of collisions. This exaggerates the effects of dis- 
sipation, producing steeper profiles. 
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